

package MGD2;

use strict;
use base qw(Exporter);

use List::Util qw(reduce);
use Math::MatrixReal;

our @EXPORT_OK = qw(
                    ComputeMGD
                );

our %EXPORT_TAGS = (
    all => [ @EXPORT_OK ],
);

use constant PI => 4 * atan2(1,1);

#
#Compute the mean value for a vector
#Parameters:
#    vector: reference to an input vector (list)
#Output:
#    the mean value for that vector
sub mean_vector {
    my $vector = shift;

    return 0 if (@$vector == 0);

    return (reduce {$a + $b} @$vector) / scalar @$vector;
}


#Compute the mean value for a matrix
#Parameters:
#    matrix: an m * n Math::Matrix,
#            m is the number of samples, n is the number of features
#Output:
#    a 1 * n one diminsion Math::Matrix
#        each value is the the mean for each feature in the matrix
sub mean_matrix {
    my $matrix = shift;

    my ($m, $n) = $matrix->dim();

    return () if ($m == 0 || $n == 0);

    my $t_matrix = ~$matrix;

    return Math::MatrixReal->new_from_rows([[map { $t_matrix->row($_)->norm_sum / $m } 1 .. $n]]);
}


#Compute the probility 
#Parameters:
#    matrix: an m * n Math::Matrix,
#            m is the number of samples, n is the number of features
#Output:
#    an m * 2 Matrix, the first column is the probility, 
#                     the second column is the offset divided by the standard deviation,
#                           if it is larger than 4, likely to be anomaly ( < 0.1% )
#
sub ComputeMGD($) {
    my $array_ref = shift;
    my $x = Math::MatrixReal->new_from_rows($array_ref);
#print "Matrix:\n", $x;

    my $mu = mean_matrix($x);
#print "Mu:\n", $mu;

    my ($m, $n) = $x->dim;

    my $sigma = Math::MatrixReal->new($n, $n); # Zero(n * n)

    for my $i (1 .. $m) {
        my $xi = $x->row($i);
        $sigma += ~($xi - $mu) * ($xi - $mu); 
    }
    $sigma = $sigma / $m;
#print "Sigma:\n", $sigma; 

    my $sigma_inv = $sigma->inverse();
#$sigma_inv->print("SigmaInv:\n");

    #the standard deviation 
    my $st_dev = sqrt($sigma->det());
#print "sd: $st_dev\n";
    my $factor = 1 / ( ((2 * PI) ** ($n / 2)) * $st_dev );

    #compute p(x) for each x
    my $probs = [];
    for my $i (1 .. $m) {
        my $xi = $x->row($i);;
        my $dist = ($xi - $mu)->norm_frobenius(); #offset from mu
        my $product = ($xi - $mu) * $sigma_inv * ~($xi - $mu);

        push @$probs, [ $factor * exp(-0.5 * $product->element(1,1)), $dist/$st_dev,
                        @{$array_ref->[$i-1]} ];
    }

    return $probs;
}


1;
